﻿#pragma once

#include <cassert>
#include <cfloat>
#include <cmath>
#include <deque>
#include <initializer_list>
#include <iostream>
#include <map>
#include <memory>
#include <stdexcept>
#include <string>
#include <type_traits>
#include <vector>

// 定义宏替换为智能指针
#define Handle(Object) std::shared_ptr<Object>

// 窗口放缩比例
// #define CAGD_SCALE_RATIO 1.25
#define CAGD_SCALE_RATIO 1.0

// 比较精度
#define CAGD_ABS_ACC 1e-5
#define CAGD_REL_ACC 1e-5

//////////////////////////////////////////////////////////////////////////////////////////
///
///

/**
 * 项目命名空间.
 */
namespace CAGD
{

//////////////////////////////////////////////////////////////////////////////////////////
///
///

/**
 * 类型枚举.
 */
enum CoordsType
{
    Vec,
    Pnt
};

/**
 * 基本数据向量.
 */
template <int N, typename T, CoordsType Type> class TCoordN
{
  protected:
    /**
     * 坐标向量.
     */
    T m_coords[N];

  public:
    /**
     * @brief 记录长度.
     *
     */
    enum
    {
        Size = N
    };

    /**
     * @brief 基本空构造.
     * @see
     */
    TCoordN()
    {
        for (int i = 0; i < N; i++)
            m_coords[i] = 0;
    }

    /**
     * @brief 基本构造（通过显式调用防止错误调用此构造）.
     * @param X
     * @see
     */
    explicit TCoordN(T X)
    {
        m_coords[N - 1] = X;
    }

    /**
     * @brief 递归构造.
     * @param t
     * @param ...args
     * @see
     */
    template <typename... Args> TCoordN(T t, Args... args) : TCoordN(args...)
    {
        static_assert(N > sizeof...(args));
        m_coords[N - sizeof...(args) - 1] = t;
    }

    /**
     * @brief 列表初始化.
     * @param l 输入列表
     * @see
     */
    TCoordN(std::initializer_list<T> l)
    {
        assert(l.size() <= N);

        int i = 0;
        for (auto it = l.begin(); it != l.end(); ++it)
            m_coords[i++] = *it;
    }

    /**
     * @brief 不同类型拷贝构造.
     * @param t
     * @see
     */
    template <int N2, typename T2, CoordsType type2> TCoordN(const TCoordN<N2, T2, type2> &t)
    {
        assert(N >= N2);
        const T2 *data = t.GetData();
        for (int i = 0; i < N2; i++)
            m_coords[i] = data[i];
        for (int i = N2; i < N; i++)
            m_coords[i] = 0;
    }

    /**
     * @brief 不同类型拷贝.
     * @param t
     * @see
     */
    template <int N2, typename T2, CoordsType Type2> void operator=(const TCoordN<N2, T2, Type2> &t2)
    {
        assert(N >= N2);
        for (int i = 0; i < N2; i++)
            m_coords[i] = t2[i];
        for (int i = N2; i < N; i++)
            m_coords[i] = 0;
    }

    /**
     * @brief 加法运算.
     * @return 相加结果
     * @param coordn 相加向量
     * @see
     */
    TCoordN operator+(const TCoordN &coordn) const
    {
        TCoordN coordn2;
        for (int i = 0; i < N; i++)
            coordn2.m_coords[i] = m_coords[i] + coordn.m_coords[i];
        return coordn2;
    }

    /**
     * @brief 减法运算.
     * @return 相减结果
     * @param coordn 相减向量
     * @see
     */
    TCoordN operator-(const TCoordN &coordn) const
    {
        TCoordN coordn2;
        for (int i = 0; i < N; i++)
            coordn2.m_coords[i] = m_coords[i] - coordn.m_coords[i];
        return coordn2;
    }

    /**
     * @brief 取反操作.
     * @return
     * @see
     */
    TCoordN operator-() const
    {
        TCoordN coordn2;
        for (int i = 0; i < N; i++)
            coordn2.m_coords[i] = -m_coords[i];
        return coordn2;
    }

    /**
     * @brief 内积运算.
     * @return 内积
     * @param coordn 相乘向量
     * @see
     */
    T operator*(const TCoordN &coordn) const
    {
        T r = 0;
        for (int i = 0; i < N; i++)
            r += m_coords[i] * coordn.m_coords[i];
        return r;
    }

    /**
     * @brief 伸缩运算.
     * @return 伸缩结果
     * @param scale 伸缩值
     * @see
     */
    TCoordN operator*(T scale) const
    {
        TCoordN coordn2;
        for (int i = 0; i < N; i++)
            coordn2.m_coords[i] = m_coords[i] * scale;
        return coordn2;
    }

    /**
     * @brief 伸缩运算.
     * @return 伸缩结果
     * @param scale 伸缩值
     * @see
     */
    TCoordN operator/(T scale) const
    {
        assert(fabs(scale) > FLT_MIN);
        TCoordN coordn2;
        for (int i = 0; i < N; i++)
            coordn2.m_coords[i] = m_coords[i] / scale;
        return coordn2;
    }

    /**
     * @brief 叉积运算.
     * @return 叉积结果
     * @param coordn 叉积向量
     * @see
     */
    TCoordN operator%(const TCoordN &coordn) const
    {
        return {m_coords[1] * coordn.Z() - m_coords[2] * coordn.Y(),
                -m_coords[0] * coordn.Z() + m_coords[2] * coordn.X(),
                m_coords[0] * coordn.Y() - m_coords[1] * coordn.X()};
    }

    /**
     * @brief 比较运算.
     * @return 是否相等
     * @param coordn 比较向量
     * @see
     */
    bool operator==(const TCoordN &coordn) const
    {
        return (coordn - *this).Modulus() < CAGD_ABS_ACC;
    }

    /**
     * @brief 自加运算.
     * @param coordn 相加向量
     * @see
     */
    void operator+=(const TCoordN &coordn)
    {
        for (int i = 0; i < N; i++)
            m_coords[i] += coordn.m_coords[i];
    }

    /**
     * @brief 自减运算.
     * @param coordn 相减向量
     * @see
     */
    void operator-=(const TCoordN &coordn)
    {
        for (int i = 0; i < N; i++)
            m_coords[i] -= coordn.m_coords[i];
    }

    /**
     * @brief 自伸缩运算.
     * @param scale 伸缩值
     * @see
     */
    void operator*=(T scale)
    {
        for (int i = 0; i < N; i++)
            m_coords[i] *= scale;
    }

    /**
     * @brief 自伸缩运算.
     * @param scale 伸缩值
     * @see
     */
    void operator/=(T scale)
    {
        assert(fabs(scale) > FLT_MIN);
        for (int i = 0; i < N; i++)
            m_coords[i] /= scale;
    }

    /**
     * @brief 自叉积运算.
     * @param coordn 叉积向量
     * @see
     */
    void operator%=(const TCoordN &coordn)
    {
        *this = *this % coordn;
    }

    /**
     * @brief 指标返回.
     * @return
     * @param i
     * @see
     */
    T operator[](int i) const
    {
        return m_coords[i];
    }

    /**
     * @brief 指标返回.
     * @return
     * @param i
     * @see
     */
    T &operator[](int i)
    {
        return m_coords[i];
    }

    /**
     * @brief 坐标 x.
     * @return coords[0]
     * @see
     */
    T X() const
    {
        return m_coords[0];
    }

    /**
     * @brief 坐标 x.
     * @return coords[0]
     * @see
     */
    void SetX(T x)
    {
        m_coords[0] = x;
    }

    /**
     * @brief 坐标 y.
     * @return coords[1]
     * @see
     */
    T Y() const
    {
        return m_coords[1];
    }

    /**
     * @brief 坐标 y.
     * @return coords[1]
     * @see
     */
    void SetY(T y)
    {
        m_coords[1] = y;
    }

    /**
     * @brief 坐标 z.
     * @return coords[2]
     * @see
     */
    T Z() const
    {
        return m_coords[2];
    }

    /**
     * @brief 坐标 z.
     * @return coords[2]
     * @see
     */
    void SetZ(T z)
    {
        m_coords[2] = z;
    }

    /**
     * @brief 坐标 w.
     * @return coords[2]
     * @see
     */
    T W() const
    {
        return m_coords[3];
    }

    /**
     * @brief 坐标 w.
     * @return coords[2]
     * @see
     */
    void SetW(T w)
    {
        m_coords[3] = w;
    }

    /**
     * @brief 取模.
     * @return 模长
     * @see
     */
    T Modulus() const
    {
        return sqrt(SquareModulus());
    }

    /**
     * @brief 模平方.
     * @return 模平方
     * @see
     */
    T SquareModulus() const
    {
        T r = 0;
        for (int i = 0; i < N; i++)
            r += m_coords[i] * m_coords[i];
        return r;
    }

    /**
     * @brief 正则化.
     * @see
     */
    void Normalize()
    {
        T mod = Modulus();
        if (mod > FLT_MIN)
        {
            for (int i = 0; i < N; i++)
                m_coords[i] /= mod;
        }
        else
        {
            for (int i = 0; i < N; i++)
                m_coords[i] = 0;
        }
    }

    /**
     * @brief 返回正则化.
     * @return 正则化向量
     * @see
     */
    TCoordN Normalized() const
    {
        T mod = Modulus();
        TCoordN coordn2;
        if (mod > FLT_MIN)
        {
            for (int i = 0; i < N; i++)
                coordn2.m_coords[i] = m_coords[i] / mod;
        }
        return coordn2;
    }

    /**
     * @brief 绝对值.
     * @see
     */
    void Absolute()
    {
        for (int i = 0; i < N; i++)
            m_coords[i] = std::abs(m_coords[i]);
    }

    /**
     * @brief 返回绝对值向量.
     * @return
     * @see
     */
    TCoordN Absoluted() const
    {
        TCoordN coordn2;
        for (int i = 0; i < N; i++)
            coordn2[i] = std::abs(m_coords[i]);
        return coordn2;
    }

    /**
     * @brief 计算夹角，返回 0-pi.
     * @return 夹角
     * @param v
     * @see
     */
    T Angle(const TCoordN &v) const
    {
        return acos(Normalized() * v.Normalized());
    }

    /**
     * @brief 考虑 this % v 与 ref 的夹角，如果夹角为正，则返回正值，否则返回负值.
     * @return 夹角
     * @param v
     * @param ref
     * @see
     */
    T AngleRef(const TCoordN &v, const TCoordN &ref) const
    {
        auto r = (*this % v).Normalized() * ref.Normalized();
        if (r > 0)
            return Angle(v);
        else
            return -Angle(v);
    }

    /**
     * @brief 到点的距离.
     * @return 距离
     * @param p
     * @see
     */
    T Distance(const TCoordN &p) const
    {
        return sqrt(SquareDistance(p));
    }

    /**
     * @brief 到点的距离平方.
     * @return 距离平方
     * @param p
     * @see
     */
    T SquareDistance(const TCoordN &p) const
    {
        return (*this - p).SquareModulus();
    }

    /**
     * @brief 在 i 分量上投影，即取消 i 分量.
     * @return
     * @param i
     * @see
     */
    TCoordN<N - 1, T, Type> Projected(int i) const
    {
        TCoordN<N - 1, T, Type> coordn2;
        for (int j = 0; j < i; j++)
            coordn2[j] = m_coords[j];
        for (int j = i + 1; j < N; j++)
            coordn2[j - 1] = m_coords[j];
        return coordn2;
    }

    /**
     * @brief 获取向量数组.
     * @return 常向量数组
     * @see
     */
    const T *GetData() const
    {
        return m_coords;
    }

    /**
     * @brief 修改向量数组.
     * @return 向量数组
     * @see
     */
    T *ChangeData()
    {
        return m_coords;
    }

    /**
     * @brief 输出函数.
     * @return
     * @param os
     * @param v
     * @see
     */
    friend std::ostream &operator<<(std::ostream &os, const TCoordN &v)
    {
        os << "{";
        for (int i = 0; i < v.Size; i++)
            os << v[i] << " ";
        os << "\b}" << std::endl;
        return os;
    }
};

//////////////////////////////////////////////////////////////////////////////////////////
///
/// 特化类型
template <int N> using BVecN = TCoordN<N, double, CoordsType::Vec>;
template <int N> using BPntN = TCoordN<N, double, CoordsType::Pnt>;
template <int N> using BVecNf = TCoordN<N, float, CoordsType::Vec>;
template <int N> using BPntNf = TCoordN<N, float, CoordsType::Pnt>;

template <class T2> using BPntT = TCoordN<3, T2, CoordsType::Pnt>;
template <class T2> using BVecT = TCoordN<3, T2, CoordsType::Vec>;
template <class T2> using BPnt2dT = TCoordN<2, T2, CoordsType::Pnt>;
template <class T2> using BVec2dT = TCoordN<2, T2, CoordsType::Vec>;

// 高精度版本
using BVec = BVecN<3>;
using BPnt = BPntN<3>;
using BVec2d = BVecN<2>;
using BPnt2d = BPntN<2>;
using BVec4d = BVecN<4>;
using BPnt4d = BPntN<4>;

// 浮点版本
using BVecf = BVecNf<3>;
using BPntf = BPntNf<3>;
using BVec2df = BVecNf<2>;
using BPnt2df = BPntNf<2>;
using BVec4df = BVecNf<4>;
using BPnt4df = BPntNf<4>;

//////////////////////////////////////////////////////////////////////////////////////////
///
///

/**
 * 模板矩阵类型.
 */
template <int M, int N, typename T> class TMatMxN
{
  protected:
    /**
     * 矩阵向量.
     */
    T m_data[M][N];

  public:
    /**
     * @brief 记录尺寸.
     *
     */
    enum
    {
        Row = M,
        Col = N
    };

    /**
     * @brief 基本空构造，默认为对角阵.
     * @see
     */
    TMatMxN()
    {
        // 对角线全一，其余归零
        for (int i = 0; i < M; i++)
        {
            for (int j = 0; j < N; j++)
            {
                if (i == j)
                    m_data[i][j] = 1;
                else
                    m_data[i][j] = 0;
            }
        }
    }

    /**
     * @brief 列表初始化.
     * @param l 输入列表
     * @see
     */
    TMatMxN(std::initializer_list<T> l)
    {
        assert(l.size() == M * N);

        auto it = l.begin();
        for (int i = 0; i < M; i++)
        {
            for (int j = 0; j < N; j++)
            {
                m_data[i][j] = *it;
                it++;
            }
        }
    }

    /**
     * @brief 加法.
     * @return
     * @param mat
     * @see
     */
    TMatMxN operator+(const TMatMxN &mat) const
    {
        TMatMxN mat2;
        for (int i = 0; i < M; i++)
            for (int j = 0; j < N; j++)
                mat2.m_data[i][j] = m_data[i][j] + mat.m_data[i][j];
        return mat2;
    }

    /**
     * @brief 减法.
     * @return
     * @param mat
     * @see
     */
    TMatMxN operator-(const TMatMxN &mat) const
    {
        TMatMxN mat2;
        for (int i = 0; i < M; i++)
            for (int j = 0; j < N; j++)
                mat2.m_data[i][j] = m_data[i][j] - mat.m_data[i][j];
        return mat2;
    }

    /**
     * @brief 取反.
     * @return
     * @param mat
     * @see
     */
    TMatMxN operator-() const
    {
        TMatMxN mat2;
        for (int i = 0; i < M; i++)
            for (int j = 0; j < N; j++)
                mat2.m_data[i][j] = -m_data[i][j];
        return mat2;
    }

    /**
     * @brief 乘向量.
     * @return
     * @param v
     * @see
     */
    BVecNf<M> operator*(const BVecNf<N> &v) const
    {
        BVecNf<M> v2;
        for (int i = 0; i < M; i++)
        {
            v2[i] = 0;
            for (int j = 0; j < N; j++)
                v2[i] += m_data[i][j] * v[j];
        }
        return v2;
    }

    /**
     * @brief 乘矩阵.
     * @return
     * @param mat
     * @see
     */
    template <int N2> TMatMxN<M, N2, T> operator*(const TMatMxN<N, N2, T> &mat) const
    {
        TMatMxN<M, N2, T> mat2;
        for (int i = 0; i < M; i++)
        {
            for (int j = 0; j < N2; j++)
            {
                mat2(i, j) = 0;
                for (int k = 0; k < N; k++)
                    mat2(i, j) += m_data[i][k] * mat(k, j);
            }
        }
        return mat2;
    }

    /**
     * @brief 伸缩运算.
     * @return 伸缩结果
     * @param scale 伸缩值
     * @see
     */
    TMatMxN operator*(T scale) const
    {
        TMatMxN mat2;
        for (int i = 0; i < M; i++)
            for (int j = 0; j < N; j++)
                mat2.m_data[i][j] = m_data[i][j] * scale;
        return mat2;
    }

    /**
     * @brief 伸缩运算.
     * @return 伸缩结果
     * @param scale 伸缩值
     * @see
     */
    TMatMxN operator/(T scale) const
    {
        assert(fabs(scale) > FLT_MIN);
        TMatMxN mat2;
        for (int i = 0; i < M; i++)
            for (int j = 0; j < N; j++)
                mat2.m_data[i][j] = m_data[i][j] / scale;
        return mat2;
    }

    /**
     * @brief 比较运算.
     * @return 是否相等
     * @param mat 比较矩阵
     * @see
     */
    bool operator==(const TMatMxN &mat) const
    {
        for (int i = 0; i < M; i++)
            for (int j = 0; j < N; j++)
                if (fabs(m_data[i][j] - mat.m_data[i][j]) > CAGD_ABS_ACC)
                    return false;
        return true;
    }

    /**
     * @brief 自加运算.
     * @param mat
     * @see
     */
    void operator+=(const TMatMxN &mat)
    {
        for (int i = 0; i < M; i++)
            for (int j = 0; j < N; j++)
                m_data[i][j] += mat.m_data[i][j];
    }

    /**
     * @brief 自减运算.
     * @param mat
     * @see
     */
    void operator-=(const TMatMxN &mat)
    {
        for (int i = 0; i < M; i++)
            for (int j = 0; j < N; j++)
                m_data[i][j] -= mat.m_data[i][j];
    }

    /**
     * @brief 自伸缩运算.
     * @param scale 伸缩值
     * @see
     */
    void operator*=(T scale)
    {
        for (int i = 0; i < M; i++)
            for (int j = 0; j < N; j++)
                m_data[i][j] *= scale;
    }

    /**
     * @brief 自伸缩运算.
     * @param scale 伸缩值
     * @see
     */
    void operator/=(T scale)
    {
        assert(fabs(scale) > FLT_MIN);
        for (int i = 0; i < M; i++)
            for (int j = 0; j < N; j++)
                m_data[i][j] /= scale;
    }

    /**
     * @brief 取值.
     * @return
     * @param i
     * @param j
     * @see
     */
    T operator()(int i, int j) const
    {
        return m_data[i][j];
    }

    /**
     * @brief 取值.
     * @return
     * @param i
     * @param j
     * @see
     */
    T &operator()(int i, int j)
    {
        return m_data[i][j];
    }

    /**
     * @brief 输出函数.
     * @return
     * @param os
     * @param m
     * @see
     */
    friend std::ostream &operator<<(std::ostream &os, const TMatMxN &m)
    {
        os << "{" << std::endl;
        for (int i = 0; i < m.Row; i++)
        {
            os << "{";
            for (int j = 0; j < m.Col; j++)
                os << "\t" << m(i, j);
            os << "\t}" << std::endl;
        }
        os << "}" << std::endl;
        return os;
    }
};

//////////////////////////////////////////////////////////////////////////////////////////
///
/// 类型特化
using BMat4x4 = TMatMxN<4, 4, float>;
using BMat4x3 = TMatMxN<4, 3, float>;
using BMat4x2 = TMatMxN<4, 2, float>;
using BMat4x1 = TMatMxN<4, 1, float>;

using BMat3x4 = TMatMxN<3, 4, float>;
using BMat3x3 = TMatMxN<3, 3, float>;
using BMat3x2 = TMatMxN<3, 2, float>;
using BMat3x1 = TMatMxN<3, 1, float>;

using BMat2x4 = TMatMxN<2, 4, float>;
using BMat2x3 = TMatMxN<2, 3, float>;
using BMat2x2 = TMatMxN<2, 2, float>;
using BMat2x1 = TMatMxN<2, 1, float>;

using BMat1x4 = TMatMxN<1, 4, float>;
using BMat1x3 = TMatMxN<1, 3, float>;
using BMat1x2 = TMatMxN<1, 2, float>;
using BMat1x1 = TMatMxN<1, 1, float>;

//////////////////////////////////////////////////////////////////////////////////////////
///
///

/**
 * 包围盒结构.
 */
struct BBox
{
    BPnt center;
    BVec x1, x2;
    BVec y1, y2;
    BVec z1, z2;
};

//////////////////////////////////////////////////////////////////////////////////////////
///
///

/**
 * 三角形结构.
 */
struct Triangle
{
    int i, j, k;
};

/**
 * 三角化结构.
 */
struct TriangleMesh
{
    std::vector<BPntf> vertices;
    std::vector<BVecf> normals;
    std::vector<BPnt2df> texcoords;
    std::vector<Triangle> triangles;
};

//////////////////////////////////////////////////////////////////////////////////////////
///
///

/**
 * @brief 组合两个点.
 * @return
 * @param p1
 * @param p2
 * @param a
 * @param b
 * @see
 */
template <class T, int N>
inline TCoordN<N, T, CoordsType::Pnt> combine(TCoordN<N, T, CoordsType::Pnt> p1, TCoordN<N, T, CoordsType::Pnt> p2, T a,
                                              T b)
{
    return (p1 * a + p2 * b) / (a + b);
}

/**
 * @brief 判断三点是否共线.
 * @return
 * @param p1
 * @param p2
 * @param p3
 * @see
 */
template <class T>
inline bool on_line(TCoordN<3, T, CoordsType::Pnt> p1, TCoordN<3, T, CoordsType::Pnt> p2,
                    TCoordN<3, T, CoordsType::Pnt> p3)
{
    TCoordN<3, T, CoordsType::Vec> vec1 = p1 - p2;
    TCoordN<3, T, CoordsType::Vec> vec2 = p1 - p3;

    // 叉积足够大，则认为不共线
    return (vec1 % vec2).Modulus() <= (FLT_MIN * vec1.Modulus() * vec2.Modulus());
}

/**
 * @brief 判断四点是否共面.
 * @return
 * @param p1
 * @param p2
 * @param p3
 * @param p4
 * @see
 */
template <class T>
inline bool on_face(TCoordN<3, T, CoordsType::Pnt> p1, TCoordN<3, T, CoordsType::Pnt> p2,
                    TCoordN<3, T, CoordsType::Pnt> p3, TCoordN<3, T, CoordsType::Pnt> p4)
{
    TCoordN<3, T, CoordsType::Vec> vec1 = p1 - p2;
    TCoordN<3, T, CoordsType::Vec> vec2 = p1 - p3;
    TCoordN<3, T, CoordsType::Vec> vec3 = p1 - p4;

    // 计算三个向量的行列式，如果足够大，则认为不共面
    return fabs(((vec1 % vec2) * vec3)) <= (FLT_MIN * vec1.Modulus() * vec2.Modulus() * vec3.Modulus());
}

/**
 * @brief 阶乘函数.
 * @return
 * @param n
 * @see
 */
inline int factorial(int n)
{
    int i = 1;
    for (int k = 0; k < n; k++)
        i *= k + 1;
    return i;
}

/**
 * @brief 排列数.
 * @return
 * @param n
 * @param i
 * @see
 */
inline int permutation(int n, int i)
{
    int j = 1;
    for (int k = n; k > n - i + 1; k--)
        j *= k;
    return j;
}

/**
 * @brief 组合数.
 * @return
 * @param n
 * @param i
 * @see
 */
inline int combine(int n, int i)
{
    return factorial(n) / (factorial(i) * factorial(n - i));
}

/**
 * @brief 计算伯恩斯坦多项式 B^n_i(t).
 * @return
 * @param n
 * @param i
 * @param t
 * @see
 */
template <class T> inline T bernstein(int n, int i, T t)
{
    T r = 1;
    for (int j = 0; j < i; j++)
        r = r * (n - j) / (j + 1);
    r = r * pow(1 - t, i) * pow(t, n - i);
    return r;
}

/**
 * @brief 笛卡尔坐标转齐次坐标.
 * @return
 * @param v
 * @see
 */
template <int N, class T, CoordsType Type>
inline TCoordN<N + 1, T, Type> cartesian_to_homogenous(const TCoordN<N, T, Type> &v)
{
    static_assert(N > 0);

    TCoordN<N + 1, T, Type> v2;
    for (int i = 0; i < N; i++)
        v2[i] = v[i];

    switch (Type)
    {
    case CAGD::Vec:
        break;
    case CAGD::Pnt:
        v2[N] = 1;
        break;
    }
    return v2;
}

/**
 * @brief 齐次坐标转笛卡尔坐标.
 * @return
 * @param v
 * @see
 */
template <int N, class T, CoordsType Type>
inline TCoordN<N - 1, T, Type> homogenous_to_cartesian(const TCoordN<N, T, Type> &v)
{
    static_assert(N > 1);

    TCoordN<N - 1, T, Type> v2;
    for (int i = 0; i < N - 1; i++)
        v2[i] = v[i];

    switch (Type)
    {
    case CAGD::Vec:
        break;
    case CAGD::Pnt:
        v2 /= v[N - 1];
        break;
    }
    return v2;
}

} // namespace CAGD
